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FLUTTER: A FINITE ELEMENT PROGRAM FOR AERODYNAMIC 

I /STABILITY ANALYSIS OF 6ENEIAL SNELLS OF REVOLUTION 
WITH THERMAL PRESTRESS 

■7 

Dennis J. Fallon 1 and Earl A. Thornton 2 
SUMMARY 

DoctiMntat ioc for tha computer progra FLUTTER it presented. A ganaral 
discussion of tha theory of aerodynamic instability with thermal prestress 
is given. Theoretical aspects of the finite elmnent matrices required in 
the aerodynamic instability analysis are also discussed. Osneral 
organisation of the computer program is explained, and instructions are then 
presented for the execution of the program. 

INTRODUCTION 

The computer progrmi FLUTTER was written in the course of research 
aimed at evaluating the effects of thermal prestressing on aerodynamic 
instability (flutter) characteristics of general shells of revolution. 
Specifically, the objective of the research was to compare the conventional 
finite element technique to the integrated finite element technique. 
Interested readers should refer to reference 1 for further details of the 
approach. 

The mein body of this report presents: (1) concept of aerodynamic 

instability, (2) finite element formulations, (3) progr as organisation and 
(4) input instruction. 


1 Assistant Professor, Department of Civil Bugineering, Old Dominion 
University, Norfolk, Virginia 23S08. 

2 Associate Professor, Department of Mechanical Engineering and Mechanics, 
Old Dominion Uhiversity, Norfolk, Virginia 23508. 


CONCEPT OP AERODYNAMIC INSTABILITY 


The free vibration equations of notion for a finite element analysis 
of shell of revolution subjected to the effect of prestressing forces 
and aerodynamic pressure is expressed as (2) : 

lK e ] {q} + tK g l {q} ♦ * [Aj {q } ♦ [M] fo} - {0 } (1) 

where the matrices [K ], [K ] , [A 1 and [M] represents the first order 

€ g G 

stiffness, initial stress (geometric), aerodynamic and mass matrices, 
respectively. The vectors {q} and {q} represent the nodal displacements 
and accelerations as a function of time. The term A represents an 
aerodynamic coefficient which is a function of the stagnation pressure and 
the Mach number. The derivation of all matrices in equation (1) is given in 
the following sections. It should be noted that prestressing of the shell 
is incorporated through the initial stress matrix. 

Now assuming that the displacement varies as harmonic function of time, 
{q} - ft} e i( * (2) 

where {q} is a vector of nodal displacement independent of time and u> is 
the natural frequency of the system, equation (1) reduces to the classical 
dynamic equation: 

( [K e ] + [K g ] + X[A e ] -w 2 (Ml) {q} e lUt - {0 } (3) 

For a solution to equation (3) to exist the determinant of the equation in 
the parenthesis oust vanish, That is 

|[K e l + [K g ] + X[A e J - ^ CM] | - 0 (4) 

The objective of an aerodynamic instability analysis is to seek a set 
of vibration inodes that are unbounded in the time domain. This criteria is 
achieved when the natural frequencies, u, defined in equation (2) are 
complex quantities. When X*0 in equation (4) the problem degenerates into 
the calculation of the in-vacuo natural frequencies of the free vibration 


case. The matrices fK ], [K ] and [M] are syamwtric, and the eigenvalues 

® 8 

are real. 

As X ia increased from aero, two of the eigenvalues approach each 
other and coalesce at a critical value of X designated A . As the value 
of X is increased beyond X the eigenvalues become complex conjugates. 

A typical plot of the natural frequencies (square root of the eigenvalues) 
versus the aerodynamic constant is illustrated in figure 1(a). Therefore, 
the value of A represents the onset of flutter of the shell. 

As will be shown later the interpolation function used in the finite 
formulation will be a Fourier series. Hence, the analysis of a shell 
reduces to seeking the eigenvalue solution of equation (4) for each har- 
monic. To obtain a complete solution to a given shell, all harmonics must 
be searched to determine the lowest value for X. This is illustrated in 
figure 1(b). Interested readers are referred to reference 3 for an expedi- 
ent technique to determine the critical harmonic. 

FINITE ELEMENT FORMULATION 
General Remarks 

The classical finite element technique involves the modeling of a lsrge 
complex system by the ass>»mblage of saialler elements (4). For the structur- 
al analyses in this report, a geometrically exact shell element was employ- 
ed. Figure 2 illustrates a typical element where Rj and R 2 are princi- 
pal radii of curvature. This element has been shown to produce excellent 
results in the computation of natural frequencies and mode shapes foi 
general shells of revolution (5). The interpolation functions for this 
element are expressed as a Fourier series in the circumferential direction 
and simple polynomials in the meridional direction as follows: 

u ■ £ TJ (s) cos n0 

n“o 

v ■ J V (s) sin n9 
n*0 n 

w ■ l W (s) cos n9 (5) 

n-0 n 
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AERODYNAMIC COEFFICIENT NATURAL FREQUENCIES 




T 



(a) AERODYNAMIC COEFFICIENT 



(b) CRITICAL HARMONIC 


Figure 1. Evaluation of critical aerodynamic coefficient. 
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FINITE ELEMENT MODEL 
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u-'IUn cosCne) 
v*I V n sinCno) 
w-I Wn cosCne) 


AXIS OF REVOLUTION 


Figure 2. Typical shell elemer.t. 
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where u, v, end w represent displacements in the meridional, circumfer- 

entiel and normal directions, respectively; 0 , V , W ere simple 

Ann 

polynomials expressed as: 


U - a ♦ a, a ♦ a, a 2 ♦ a. s 

n l,n 2,n 3,n 4,n 


V ■ a, ♦ a, s + a. s^ ♦ a. s^ 
n 5,n 6,n 7,n 8,n 


W 

n 


( 6 ) 


where a. represents the ith generalised coordinate for the nth harmonic; 

i ,n 

s is the meridional coordinate and 0 is the circumferential coordinate. 

At each nodal circle designated by 1 and j in figure 2 , there are 

seven degrees of freedom: w, u, v, (w' - “•), u' , v' , (w*‘ - ■■[—). The 

*1 ” 1 

prime denotes the differentiation with respect to s. 


Pirst Order Strain Energy 

The first order stiffness matrix defined in equation (1) can be derived 
from the first variation of the first order potential energy. The potential 
energy for a general shell of revolution is: 

1 2 22 
V * — JJ (C £j ^ 2C |2 £2 ^ ^22 ^2 ^ Cgg £ d6ds 

2 

1 ..2 2 2 
+ ~ JJ '^11 + *1 *2 + ^22 *2 + ^66 *12^ r d0ds (7) 

2 

where £j, e 2 , £3, are first order meridional, circumferential and 

shearing strains according to Novoshilov shell theory (6), Kp k 2 , tc 12 
are first order meridional, circumferential and cross curvature according to 
Novoshilov shell theory; C fcl and D fcl are elastic stiffness coefficients. 
By use of the definitions (see Appendix A) of the strains and curvature 
which are explicit functions of the displacements and, with the additional 
use of the orthogonality properties of the sine and cosine, i.e.: 


6 


0 


a*n; «Fn*0 
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2* 

/ tin i)6 tin a6d9 » 
o 


w 


■FD 


2» 

/ cot n6 cot «6d6 * 
o 


0 

2w 

« 


a<0 

m»n«0 

a"n*0 


(B) 


the first order stiffness sitrix in tsrns of the nth harmonic is obtained. 
Exset expression for this aatrix esn be found in reference S. Note thet 
integrstion with respect to the meridional coord inete s wee performed 
using s ten point Gsuss achsae. 


Second Order Strsin Energy 

A consistent initisl stress (geoaetric) astrix which incorporates the 
effects of prestressing is foraulsted froa the contribution of the socond 
order strains in the strain energy expression. This energy is (7): 


U 

( 


N* // «j 2) r d6ds ♦ N* e // t[ 2) r d6ds 


♦ 12 


(w ee - * *:, > 

<l-v*> 


II U*1 V « 2 J V f l ~ v l *8 C i2^ ” 


d fids 


(9) 


t* co» ♦ 


where are initial stresses and N* # , Mgg are initial ao a ents 

due to prestressing in the swridional and circuaferential directions, 

(2) (2) 

respectively; Cj , c 2 are second-border strains, 6 is the slope of the 
shell surface in the meridional direction; and Bg, 6e * r * perturbation 
rotations in the meridional and circuaferential directions, respectively. 

The developaent of the initial stress aatrix follows in a siailar manner as 
the first order stiffness aatrix. For definitions of the strains and 
rotations in taras of the displaceaents refer to Appendix A. Integration 
with respect to the aeridional coordinate was also performed by a ten point 
Gauss scheaa. 


Kinetic Energy 

The consistent mass matrix used in this study was derived free the 
kinetic energy of the system. Specifically, the kinetic energy is expressed 
as (5) : 
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E . £ // {'2 «. y 2 + ; 2 } t r dfldi ( 10 ) 

k 2 

where p ia the seas denaity; t ia the ahell thickneaa end e dot denotea 
differentiation with respect to time. 

Aerodynamic Virtual Work 

The aerodynamic matrix expressed in equation (1) was derived from the 
virtual work of the aerodynamic forces acting on the ahell. These aero- 
dynamic forces are formulated using a first order, high Mach number (>1.7) 
approximation to linear potential flow theory ( 8 ): 

«W - - // 6 wX [— ♦ - l' M “— ) v - 3 ^] r d Ms (11) 

3s U - 

7a „ 1/2 

where \ ■ 8 « (M‘ - 1) ; is the freestream Mach number, q is the 

freestream dynamic press re and U ic the local flow velocity, For exact 
expression for the aerodynamic matrix refer to Appendix B. 

PROGRAM ORGANIZATION 


The computer program consists of four basic steps: (1) reading of 

input data, ( 2 ) the evaluation of element matrices, (3) the arsembly of 
element matrix to form g/.obal matrices and (4) the evaluation of eigenvalues 
and eigenvectors. A flow chart of *:ne program is illustrated in figure 
3. 

The input data which will be defined in more detail later in this 
report consists of material properties and geometry data of the particular 
shell. Such material properties as the modulus of elasticity and Poisson's 
ratio are input parameters for the stiffness calculation. The mass density 
is required for the evaluation of the mass matrix. These parameters are 
assumed to be constant throughout the shell. Some geometry data required on 
input is the length of each element and thickness of the shell (also assumed 
to be constant). At present the program is set up to do conical shells. To 
evaluate the flutter boundaries of other shells the subroutine Radius must 
be revised. Interested readers are referred to reference 5 for a detailed 
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Figure 3. Flow chart FLUTTER 
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procedure on how chit aty bo dono. Alto on input parameter for tho boundary 
condiciont mutt bo defined. Thit it done through the boundary condition 
code ICASE which hat the interpretation ahown in table 1. A free boundary 
condition meant all diaplacementa and rotationa are releaeed. A freely 
aupported condition impliea that ou^y the normal diaplacementa and the rota- 
tion are reletted, whereat in the clamped condition the thell ie aeaumed to 
be completely fixed. 

The calculation of the member element begint with the evaluation of the 
trana format ion matrix per each element. Thit matrix definea the relation- 
ahip between the genertliaed coordinate (tee equation and the diaplace- 
menta and rotationa at the ^nda of each element Then per each htrawnic the 
maaa, firat order atiffneaa, aerodynamic and initial atreaa matricee are 
calculated per element. Internal atreaaea are calculated at each integra- 
tion point by a thermal atreaa computer program. Theae atreaaea are read 
off Tape 11 by FLUTTER. 

After element matricea are formed they are auperimpoeed to form the 
global matricea. Rowa and columna are eliminated according to the applica- 
ble boundary condition apecified by ICASE. Then uaing a NASA/Langley 
Computer Center subroutine complex eigenvaluea and eigenvectora are computed 
per each apecified value of aerodynamic coefficient. Thia procedure la 
repeated per each hanaonic apecified. 

INPUT PROCEDURE 

The following ia the procedure per line of input data for the proper 
execution of the program. Note all input ia free formatted. 

Firat Line: IDEN 

XDEN " any alphanumeric charactera to define a particular run. 
Second Line: K, NBEC, KLAST, ICASE 

K ■ nuaber of element ; 

*BEG ■ beginning harmonic; 

NLAST " atopping harmonic; 

ICASE * boundary condition ^ode. 

Third Line: SO 

SO ■ origin of ahell'a coordinate ayatem 
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Table 1. Boundary Condiciona Codes 


Value I CASE 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 


Type of Boundary Condifcon 

Free - Free 

Free - Freely Supported 
Freely Supported - Free 
Free ~ i> imply Supported 
Simply Supported - Free 
Free ~ Clamped 
Clamped - Free 

Freely Supported - Freely Supported 
Simply Supported - Simply Supported 
Clamped - Clamped 

Freely Supported - Simply Supported 
Freely Supported - Clamped 
Simply Supported - Freely Supported 
Simply Supported - C lamped 
Clamped - Freely Supported 
Clamped - Simply Supported 


Fourth Line: E(I) L - l, K 

E(I) * length of the ith element; K values required on this 
line 

Fifth Line: Y0UNG1, Y0DNG2, XMU1, XMU2, RHO, TH, G12 

YOUNGl ■ Young's Modulus in meridional direction; 

Y0UNG2 ■ Young's Modulus in circumferential direction; 

XMU1 * Poisson's ratio in meridional direction; 

XMU2 « Poisson's ratio in circumferential direction; 

RHO “ mass dens .cy of shell; 

TH • thickness of shell; 

G12 " shear modulus of shell 

Sixth Line: MACH, NLAMB, IBUCK 

MACH ■ reference Mach number; 

NLAMB * number of aerodynamic coefficients; 

IBUCK ■ 0 (aerodynamic stability problem) 

■ 1 (bifurcation buckling loads - eigenvalues 
represent buckling loads) . 


Seventh Line : LAMB 

LAMB “ aerodynamic coefficient (This line will be repeated 
NLAMB times). 


CONCLUDING REMARKS 


A finite element prograa for the computation of the aerodynamic insta- 
bility of general shells of revolution with thermal prestress is described. 
The theoretical formulation of the finite element matrices is discussed and 
input instructions for execution of the program are given. Application of 
the program are presented in reference 1. 
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APPENDIX A: EQUATION FOE STRAINS AND NOTATION 


Definition of itraioi and rotations according to Novoshilov shell 
theory. 


£ i-»' *T7 


TT ♦ u ♦ — 

2 r 30 r R, 


e i2 • 7 13 * »’ ■ 7* ’ 


l R i 

C 1 ■ ' *" *TT --7 


■ m — . 

1_ 

S 2 * 4 1 3v _ r 

1 w* r* 

2 

r 2 

302 * R 2 80 ’ 

r r Rj 


. J, 

j 2 * ♦ i_ r t + 

1 3u + •, 

12 “ ‘ 

r 

3s 30 r 2 36 

r R x 36 1 

< 2 > - 

1 

*■ R 1 8 r 2 

[ i“ - vr* - 
l 30 

1 

2 

< 2) - 

1 

_ ri* - "f ♦ _ 

- Lsa d J 

— fcs ~ vr' 


j r2 . 

L!* - “_1 

l 3s R^ 

_1 rjhr vr -I 
r l 38 “ R 2 J 




(') denotes different at ion with respect to s. 
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APPENDIX B: AERODYNAMIC COEFFICIENTS 


Aerodynamic coefficients in terms of generalised coordinates. 



